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Motivated by the increasing interest in models which consider scalar fields as viable dark matter 
candidates, we have constructed a generalization of relativistic Boson Stars (BS) composed of two 
coexisting states of the scalar field, the ground state and the first excited state. We have studied 
the dynamical evolution of these Multi-state Boson Stars (MSBS) under radial perturbations, using 
numerical techniques. We show that stable MSBS can be constructed, when the number of particles 
in the first excited state, A'^^^', is smaller than the number of particles in the ground state, A^'^^ 
On the other hand, when A'^'^' > A'''^', the configurations are initially unstable. However, they 
evolve and settle down into stable configurations. In the stabilization process, the initially ground 
state is excited and ends in a first excited state, whereas the initially first excited state ends in a 
ground state. During this process, both states emit scalar field radiation, decreasing their number 
of particles. This behavior shows that even though BS in the first excited state are intrinsically 
unstable under finite perturbations, the configuration resulting from the combination of this state 
with the ground state produces stable objects. Finally we show in a qualitative way, that stable 
MSBS could be realistic models of dark matter galactic halos, as they produce rotation curves that 
are flatter at large radii than the rotation curves produced by BS with only one state. 



PACS numbers: 04.40.-b,04.40.Dg,95.35.-|-d 



I. INTRODUCTION 

The existence of Dark Matter (DM) in the Universe 
is strongly supported by astronomical observations that 
range from galactic up to cosmological scales (see for ex- 
ample [H and references therein). Observations indicate 
that stars rotate too fast around the center of the galaxy 
to be bound by Newtonian gravity if all matter is vis- 
ible i i. This issue, knovifn as the Rotation Curves 
(RC) problem, implies vi^ithin the context of Einstein's 
General Relativity, that a great amount of the matter in 
the galaxy is invisible. The nature of this dark matter, 
that has a negligible interaction with the visible matter 
and whose presence is only observed trough its gravita- 
tional effects, is still unknown. The most popular candi- 
date are the so-called weakly interacting massive particles 
(WIMPs) leading to the standard Cold Dark Mat- 

ter model. This scenario is very successful at a cosmolog- 
ical level, as its predictions are in good agreement with 
the observational data [1,I3- However, it has difficulties 
in fitting the observations at a galactic level [sl-fl^. If 
DM is modeled by WIMPs, one obtains a cuspy density 
profile of the DM in the galaxy. But high resolution data 
of low surface brightness galaxies, which are composed 
mainly of DM, imply that their DM distribution has a 
flat core [13, El • This model also fails in predicting the 
number of satellite galaxies around each galactic halo, 
exceeding far beyond what is observed around the Milky 
Way [11. 

A different approach consists in describing the dark 
matter as a scalar field [l3l-[l5i. The Scalar Field Dark 



Matter (SFDM) model has been proved to be successful 
at cosmological scales [l6l| . This model can also avoid the 
problems that WIMPs present at a galactic level, produc- 
ing a non-cuspy density profile [imlq and explaining the 
dearth of satellite galaxies around each galactic halo [3 ■ 
Because of the viability presented by the SFDM model, 
it is stimulating to go further on testing it. For instance, 
the model has to reproduce the observed RC of galax- 
ies. At this point Boson-Star-like objects could play an 
important role. 

In the SFDM model, the dark matter particle is 
an ultra-light massive spinless boson (m ~ 10^^'^eV 
[T6|). These bosons could collapse forming gravitation- 
ally bounded structures. With such ultra-light mass, the 
boson's Compton wave length is of the order of kilo- 
parsecs, and structures with comparable length scales - 
like galactic halos - could be formed as condensates de- 
scribed by a coherent scalar field [19]. These condensates 
can be associated with Boson stars (BS), which are solu- 
tions of the Einstein-Klein- Gordon equations where the 
gravity attraction is balanced by the dispersive charac- 
ter of the scalar field [l^, BS were first studied by 
Kaup [ll and a year later by Ruffini and Bonazzola [23| , 
who settled two different treatments. The first one, de- 
veloped in , is a completely classical treatment with a 
massive complex scalar field minimally coupled to grav- 
ity. In the second one [11], a real quantized scalar field 
is introduced in order to describe a many boson system 
though maintaining the geometry as a classical entity 
(i.e., a semiclassical limit is adopted). The relevant quan- 
tity computed in this case is the mean value (QlT^^lQ) of 
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the energy momentum tensor operator, where \Q) is the 
state of the system of many particles. When one consid- 
ers a \Q) for which all the particles are in the same state, 
it turns out that the mean value of T^i^ generates the 
same energy-momentum tensor as the complex classical 
field and consequently, the same macroscopic results, i.e. 
both a quantized real scalar field and a classical complex 
scalar field yield to the same self-gravitating system. 

Until very recently, only boson stars with all the par- 
ticles in one state have been considered. However, New- 
tonian configurations with scalar fields coexisting in the 
ground and excited states were introduced in [23| in order 
to model dark matter halos. Previous studies in model- 
ing dark matter halos using newtonian BS were done in 
However, these structures can not account for a 
realistic halo since the configurations in the ground state 
produce RC which are not flat enough at large radii. The 
case of a massless scalar field used in order to fit rota- 
tion curve data of several galaxies was considered in [sO] . 
Nevertheless, it has been proved in [3l|, HI] that no non- 
singular self-gravitating solitonic objects can be formed 
with a massless scalar field. On the other hand, RC from 
excited BS are in better agreement with the astrophysical 
observations, but these structures are unstable [33j. 

In it was shown that Newtonian mixed configura- 
tions could account for more realistic DM halos, as they 
are stable and could fit with better agreement the ob- 
served RCs even at large radii. In the present work, 
we are interested in the fully relativistic generalization 
of these Multi State Boson Stars (MSBS). The idea is 
to consider the possibility that the bosons are not all in 
the same state, but rather populating different coexisting 
states, as was already pointed out in 2i\. It turns out 
(see [23] and Appendix |A| that the resulting equations 
for a MSBS in the semiclassical approach are equivalent 
to the case where a collection of complex classical scalar 
fields is considered, one for each state, which are only 
coupled through gravity. Without any loss of generality, 
we can choose either the semiclassical formulation of such 
MSBS or its pure classical counterpart. We will follow 
the latter approach in order to investigate MSBS, analyz- 
ing in detail the properties and stability of configurations 
with two states, a ground state and a first excited state. 

The previous stability studies of BS can be divided 
roughly in two categories, depending on the type of per- 
turbations considered: 

1. Studies where the perturbations preserve the num- 
ber of particles (infinitesimal perturbations) , which 
generally involve linear perturbation analysis [sH . 
[bgI ] and catastrophe theory jSTj . 

2. Those where the perturbations don't conserve the 
number of particles (finite perturbations), which 
have been addressed mainly by numerical studies, 
[s^m,!!^. Furthermore, the late time evolution of 
unstable boson star under finite perturbation can 
only be followed by numerical simulations. 



A consistent result coming from both type of studies is 
that BS in the ground state are stable against pertur- 
bations if the amplitude of the scalar field at the origin 
0(0) is smaller than the critical value 4>max{Q) where the 
maximum mass Mmax is reached. In the case of excited 
BS there are some important differences. Although the 
linear stability analysis shows stability up to the critical 
value (pmaxiO) 21, 36] when the number of particles is 
not conserved, excited BS are intrinsically unstable even 
for (j){0) < (j)max{0), since finite perturbations drive the 
star either to collapse to a black hole or to decay to the 
ground state [H, [Sg] . 

From these results one could infer that the MSBS 
states would be unstable under perturbations when the 
number of particles is allowed to change, since they con- 
tain at least one excited state. Quite surprisingly, our 
numerical analysis shows that there is a region of the solu- 
tion space with stable configurations. Roughly speaking, 
the ground state produces a deeper gravitational poten- 
tial which can be enough to stabilize the excited state. 

This paper is organized as follows. In Section |lll we 
present the formalism used for the numerical evolution 
of the Einstein-Klein- Gordon system, describing a MSBS 
in the classical approximation. Section [IIII describes how 
the initial data for a MSBS with two different states is 
constructed. In Section IIV A| we present numerical re- 
sults obtained from the evolution of the two-state boson 
stars. We study two features of the evolution, namely 
the stability and the late time behavior. In both cases, 
we add an small perturbation. In the first case, the per- 
turbed MSBS is evolved for short time scales, in order to 
study the behavior of the perturbations. In the second 
case, we evolve unstable MSBS for longer times follow- 
ing the properties of the resulting configurations. We 
compute in section |V] RC from stable MSBS and dis- 
cuss qualitatively why these RC are in better agreement 
with the observed RC of galaxies. We conclude in Sec- 
tion |Vll The description of the semiclassical approach is 
presented in Appendix \^ while Appendix [B] is devoted 
to a detailed description of the evolution equations used 
for numerical evolution. 



II. THE EINSTEIN-KLEIN-GORDON SYSTEM 

Let us consider a semi-classical real massive scalar field 
with P different excited states, which is equivalent to con- 
sidering a collection of P classical complex scalar fields 
(one for each state) coupled only through gravity. In a 
curved spacetime, the dynamics of these MSBS can be 
described by the following Lagrangian density (adopting 
geometrical units, i.e. G = c = h = 1), 



(1) 



where R is the Ricci scalar, Qah is the spacetime met- 
ric, (Z)^"-* are the scalar fields, (^^"^ their complex con- 
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jugate, and T^d^*-"-*!^) a potential depending only on 
|^(n)|2_ Throughout this paper, Roman letters from the 
beginning of the alphabet a, 6, c, .. denote spacetime in- 
dices ranging from to 3, while letters near the middle 
I, j, k, .. range from 1 to 3, denoting spatial indices. This 
Lagrangian gives rise to the equations determining the 
evolution of the metric (Einstein equations) and those 
governing the scalar fields behavior (Klein-Gordon equa- 
tions). 

The variation of the action associated with the La- 
grangian ([T]) with respect to the metric gab, leads to the 
well-known Einstein equations 



R 



gab — Si-KTab, 



(2) 



where Rab is the Ricci tensor. Tab is the total stress- 
energy tensor, given by the addition of the single stress- 
energy tensors of each scalar field, namely 



Tab = 



^ ab 
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The Einstein equations form a system of 10 non-linear 
partial differential equations for the spacetime metric 
components gab- 
On the other hand, the variation of the Lagrangian 
dl]) with respect to each scalar field leads to a set 
of Klein-Gordon (KG) equations which are only coupled 
through the gravity, 

where the box □ = g°'''S/a^b stands for the wave operator 
on a curved background. In the following, we will restrict 
ourselves to the free field case, where the potential takes 
the form 



(6) 



with m a parameter that can be identified with the bare 
mass of the field theory. 

The matter Lagrangian is invariant under global U(l) 
transformations 



0(n) ^ ^(»)e*'^*' 



(7) 



This symmetry implies that there is a set of Noether 
currents densities Ja 



(n) 



/(«) - 



(8) 



(n) _ 



0. 



satisfying for each n the conservation law V°Jq 
The Noether charge contained in some radius is given by 



5"" 4") dx^ 



(9) 



so that the total Noether charge of the system N is the 
sum of the total individual ones iV(") = Af(")(oo), namely 



(10) 



As discussed in [23|, this quantity N can be associated 
with the total number of bosonic particles. Consequently, 
can be interpreted as the number of particles in the 
state labeled by n. 



III. INITIAL DATA FOR MULTI-STATE BOSON 
STARS 

The initial data for the MSBS configurations is com- 
puted in spherical symmetry with a one-dimensional 
code. We adopt the following harmonic ansatz for each 
scalar field. 



(11) 



With this assumption, the source for the Einstein equa- 
tions becomes time independent. Our goal is to find 
{(/)„(r), Wn} and the metric coefficients, such that the 
spacetime generated by this matter configuration is 
static. 

We begin by considering the problem in polar-areal co- 
ordinates [2^ . I4Q1] . The line element in these coordinates 
takes the form 



-a [rf dt^ + a {rf dr^ + r^dQ^ 



(12) 



Then the equilibrium equations, obtained by substituting 
the ansatz ([TT|) and the metric ((T^) in the Einstcin-Klein- 
Gordon system (|2l5p . are given by 



drCL 



drCX = — 
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,(13) 
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r 

p 

E 

n=l 



2 1 2 j,2 

m a 0„ 



= -<1 




m ] (j)na . 



(16) 



In order to obtain a solution of this system, we pro- 
vide the following boundary conditions, motivated by the 
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physical situation under study, 



0n (0) 




(0) 


= 0, 


a(0) 


= 1, 


lim (t>n (r) 


w 0, 


r— ^cxD 


lim a (r) 

r— >oo 


= lim 



a[r) ' 



(17) 
(18) 
(19) 
(20) 

(21) 



which guarantee regularity at the origin and asymptotic 
flatness. For given central values of the fields {0cn}, we 
only need to adjust the eigenvalues {a;„} and the value 
q;(0) in order to generate a solution with the appropriate 
asymptotic behavior (|20j|21[) . This is a shooting prob- 
lem that we solve by integrating from r = towards the 
outer boundary r — rout, with a second order shooting 
method. The boundary conditions for the scalar fields 
at Tout are imposed considering that localized solutions 

decrease asymptotically as 4>n ^ exp (^—^Jmi? — uj^'t^ jr 

in a Schwarzschild-type asymptotic background. At the 
outer boundary, the conditions are 



0n {^c 



1 



0. (22) 



The shooting procedure is performed for different values 
of Tout- As Tout is increasing, the shooting parameters 
converge, and we choose the solution as the one which 
satisfies the conditions ^F2\ for some Tout within a pre- 
scribed tolerance. From this point on, we match to the 
scalar fields and the metric coefficients their asymptotic 
behavior. 

A qualitative characteristic of the radial functions 0„ 
is their number of nodes (ie, how many times they do 
cross zero), which reflects the excited state of the boson 
star. If the radial function does not have any node, the 
boson star is in the ground state. When there is a node, 
the boson star is in the first excited state, and so on. 
In the next subsection we construct initial configurations 
with two scalar fields P — 2, one in the ground state and 
the other in the first excited state. Notice that this is the 
simplest non-trivial configuration, since the MSBS with 
two scalar fields in the ground state can be reduced to one 
scalar field solution by redefining the scalar fields. This 
is a consequence of the indistinguishably of the boson 
particles in the same state. 

Once the solution is computed in this coordinate sys- 
tem, a change of coordinates is performed to maximal 
isotropic ones. 



(r) dt^ + V/ (r) {dr'' + f^dil^) , (23) 



which are more convenient for our numerical evolutions. 
Finally, a simple inspection on the system (|13m6p shows 
that the re-definition 



rm , ujr, 



m 



(24) 



leads to a set of equations that no longer have m on it. 
Hence, the selection of geometrical units (e.g. G = H = 
c = 1) and the redefinition ([M)) give us dimensionless 
units for f,w„, which is equivalent to choosing m = 1 
in our equations. Throughout this paper, we use these 
dimensionaless coordinates. 



A. Configurations of ground and excited states 
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FIG. 1: Ground-lst excited configuration for 0i(O) = 0.0197 
and fraction r] — 1. The upper panel corresponds to the initial 
profiles of the two scalar fields, and the lower panel, to the 
lapse function a and the conformal factor >1/. 

Let us consider the simplest non-trivial case with only 
two scalar fields P = 2, one with TV^^^ particles in the 
ground state, and the other with N^^^^ particles in the 
first excited state. A useful way to construct the initial 
data is specifying the fraction between the number of 
particles in each state of the configuration. 



V = 



(25) 



In this case, we complete the system 



"T6| with the 

differential expressions for the number of particles in each 
state 



5^7V(")(r 



2„2 



47r-a;„0„r 
a 



(26) 
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with boundary conditions, 7V^"^(0) = 0. If 77 is speci- 
fied, it is sufficient to prescribe as boundary conditions 
the central value of only one of the scalar fields, for 
instance, (jfci- The new system of equations (j22p and 
(|26p becomes a shooting problem for the four parame- 
ters {cJi,a;2,a(0),(/)2(0)}. For a specific fraction i], it is 
necessary to adjust the four parameters such that equa- 
tion (gH) and the condition A^(^H^o«i) = V^^^K'^out) are 
satisfied. 

Fig. [T] shows an example of the radial profiles of the 
two scalar fields, lapse and conformal factor, for a MSBS 
with {tj = 1,M0) = 0.0197}. 

Two important characteristics of MSBS are the total 
gravitational mass M and the radius Rgg. The first one 
is calculated as 



M 



rout A _ 1 
2 V a?{rout)) 



(27) 



and the radius i?g9 is defined as the radius where M 
reaches the 99% of its value. The choice of geometrical 
units, plus the re-definition of the coordinates eq. (1^^ . 
imply that both iJgg and M are dimensionless variables. 
The physical units can be recovered by using the follow- 
ing relations: 



m 



R, 



Physical 
99 



R 



•99- 



(28) 



where mp the Planck's mass and m the mass of the boson 
associated to the scalar field. 

In Fig. [2] we have plotted the values of these two quan- 
tities for all the constructed MSBS initial configurations. 
On the top panel they are shown as functions of the cen- 
tral value of the scalar field in the ground state (0) and 
the Noether fraction 77, while in the bottom panel they 
are plotted as functions of the central value of the scalars 
fields. 

These figures already show some differences between 
single BS and MSBS configurations: there are an infinite 
number of possible equilibrium configurations (i.e. so- 
lutions for the static EKG system (1131116^ ) between the 
two extreme cases N'^'^^ = and N'-^^ = that corre- 
spond precisely to the ground state BS and first excited 
state BS respectively. Fig. [5] shows a more complex 
behavior of the mass M and Rgg than for a single BS, 
which can be observed easily in Fig. [3] where some slices 
01 (0) = constant and 02(0) = constant of these surfaces 
are displayed. After constructing these configurations, 
the next problem is stability. It is clear that the known 
results of BS stability are not immediately applicable to 
MSBS. 

The stability of the multi-states configurations under 
finite perturbations will be addressed in the next section. 



IV. NUMERICAL SIMULATIONS 

In this section, we present a numerical analysis of the 
dynamical properties of MSBS, focusing on long-term 
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FIG. 2: Total gravitational mass {M{(f>i,r]) ,Rgg{(j>i,ri), 
M{<j)i{0),(j>2{0)) and Rgg{(j>i{0), (j>2{0))) for initial data of 
MSBS configurations witli two states, the ground and the 1st 
excited state. 
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FIG. 3: Difi'erent slices of the Mass (upper panel) and the 
radius Rqq (lower panel) as a function of the central value 
of the scalar field. On the left there are slices at 02 (0) = 
constant while on the right there are 0i(O) = constant. 



stability and the final state of the unstable configura- 
tions. To this purpose, we write the Einstein-Klein- 
Gordon system as a set of evolution equations for the 
scalar fields and the metric components. We consider 
a generic spherically symmetric spacetime with the line 
element 



ds-" 



-a dt + grrdr + r gggd^l 



(29) 



where a is the lapse function and {grr, gee} are the metric 
components. Notice the explicit dependence on the factor 
r^, such that the component ggg is regular at the origin. 
This is a necessary condition in our implementation in 
order to deal with the coordinate singularity at r = 0. 

The evolution equations for the geometry are obtained 
by substituting the metric coefficients ([^^ into a partic- 
ular formulation of the Einstein equations. In this study 
we have considered the Z3 formulation [41], which in- 
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eludes the momentum constraint into the evohition sys- 
tem, by considering an additional vector Zi as an evolved 
field. Further details regarding the Z3 system in spher- 
ically symmetry can be found in [4^ . while the regular- 
ization of the coordinate singularity r = is similar to 
the one described in [43|. 

In spherical symmetry, there are independent evolu- 
tion equations only for the lapse a, the metric compo- 
nents {5rri3ee}j the extrinsic curvature {if^ji^Tg} and 
the Z- vector components \^Z^,Zq\. All these evolution 
equations are prescribed by the Einstein equations, ex- 
cept the one corresponding to the lapse, which is related 
to the choice of coordinates and can be specified freely. 
A common choice, which give rise to a hyperbolic system 
of equations, is the harmonic slicing 



dta — ~o?trK, 



(30) 



where trK = K; + 2 K^. 

On the other hand, the evolution equations for the 
scalar fields are obtained by substituting the spherically 
symmetric metric (1^^ in the Klein-Gordon equations ([5]). 

A first order reduction in space can be performed by 
introducing as independent quantities the spatial deriva- 
tives of the metric and new fields related to the time and 
spatial derivatives of the scalar field, namely 



—Ora 
a 



rr 

r y o 



- d p, 

= —Orgee, 



-dt 



a 



(31) 



In this way, we obtain a fully first order system of 
evolution equations for the geometry and the scalar fields, 
with the following set of evolution variables 

{a, grr, gee, K, kI A„ B„\ D^e' , Zr, 0^") , '/'i"^}- 

This first order system can be written in balance law form 

dtU + dk''F{V)=S{lJ), (32) 

which allows the use of advanced numerical methods 
based on Finite Volume algorithms. Details about the 
exact form of the evolution equations can be found in 
Appendix [B] 

We have implemented the equations using the Method 
of Lines, in order to separate the time and the spatial 
discretization. The time integration is performed with 
a third order Strong Stability Preserving Runge Kutta 
method [l^]. The spatial discretization is based on a 
standard fourth order centered finite difference scheme, 
plus third order accurate dissipation j^ . 



A. Stability of the MSBS 

The stability of MSBS configurations is a basic require- 
ment for considering them suitable models of galaxy ha- 
los. For a single boson star, the stability has been previ- 
ously studied both analytically and numerically, showing 



that it is stable if 0(0) < (pmax (0). In this section, we ana- 
lyze the stability of MSBS in the range (piiO) < (f'maxiQ), 
which are stable for ?; = 0, avoiding this way the too- 
massive unstable MSBS. The stability will be studied 
following only a numerical approach, by perturbing the 
MSBS and studying the evolution of this perturbation. 
We will restrict ourselves to study stability against spher- 
ically symmetric perturbations by using the equations 
described in the previous subsection. Notice that this 
is only a necessary condition for the most general case, 
since asymmetric perturbations may still be unstable. 

In order to study numerically the stability of the MSBS 
configurations, we perform the following steps: 

• Construct different initial data sets for MSBS with 
a given 0i(O), by varying the Noether fraction 77. 

• Add a real scalar field far outside the radius Rgg of 
the MSBS, which will be coupled to the MSBS only 
through gravity. The energy density corresponding 
to the scalar field is only 0.01% of the total energy 
density of the MSBS, so it will act just as a small 
perturbation with negligible errors in the form of 
constraint violations. 

• Perform evolutions of the Einstein-Klein-Gordon 
system and study the behavior of the MSBS. The 
scalar field perturbation will fall into the MSBS and 
later disperse to infinity. The gravitational inter- 
action during that time is expected to excite the 
unstable modes, if any. In this way, the modes are 
excited sooner than only by numerical errors. 

• Bracket the MSBS which lead to significant ex- 
ponential growing modes. For the stable MSBS 
configurations, the perturbations will only oscillate 
without growing. We expect that the MSBS with 
low rj will be stable, since the major contribution to 
the complete configuration comes from the stable 
ground state, while those with high 77 will corre- 
spond to unstable MSBS. 

• Fit the growth rate of the unstable MSBS for each 
set of stars with the same 0i(O) by varying 77. Ex- 
trapolate to find the maximum allowed Noether 
fraction rjmax which separates the stable and un- 
stable states. This procedure allows us to obtain 
reliable estimations, without evolving every config- 
uration in order to obtain rjmax- Moreover, it might 
be difficult to distinguish stable from unstable con- 
figurations when they are close to rjmax, since the 
exponential growth is very low in that region. 

We are going to restrict the numerical stability 
analysis to only three different values of 0i(O) — 
{0.0143,0.0197,0.0423}. Fig. |2 displays the total mass 
and the radius Rgg for these configurations, as a function 
of 77. In the simulations with 77 < 1.2, we did not detect 
any unstable exponentially growing mode, or they were 
difficult to measure for some families of solutions. The 
results indicate as the upper bound rjmax < 1.2. 
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FIG. 4: The Mass (upper panel) and the Radius 7?9g (lower 
panel) are presented as a function of -q, for three different 
central values of the ground state field 1^1(0). The points are 
the configurations used for the fitting of rjmax- 



In fig. \5\ we show the typical behavior for one of these 
simulations, corresponding to stable and unstable MSBS. 
The perturbation has an exponentially growing behavior 
only for the unstable MSBS. The maximum of the scalar 
field (j)2{r — 0) can be fitted with a function of the type 



A exp((7t) cos{ujt + ip), 



(33) 



which allows us to compute the exponential growth rate 
a. 

We performed fits for the unstable MSBS perturba- 
tions with T] > 1.2, marked with filled geometrical shapes 
in fig. m The results for a are represented in fig. [71 with 
the extrapolation to the rjmax, which in principle could 
be a function of 4>i{0). The three different families of 
solutions point to rjmax ~ 1- 

In order to show the robustness of these results we 
have repeated the simulations with more resolution and 
with a different amplitude of the perturbation. The re- 
sults are almost identical, as shown in figlHl The only 
significant difference is the unstable case with larger per- 
turbation amplitude. The perturbation seem to interact 
non-linearly with the star, and the unstable exponentially 
growing mode is excited later. In spite of this delay, the 
growth rate is identical to the other cases with smaller 
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FIG. 5: Maximum of the central value of the scalar field 
in the excited state </>2, for two different values of rj with 
4>i{0) = 0.0197. The MSBS with 7] = 0.4 is in the sta- 
ble branch and the induced perturbations do not grow. The 
MSBS with = 1.6 is clearly unstable and the perturbations 
exhibit an exponential growth. 



perturbation and higher resolution. 



B. Fate of the unstable states 

Another question which arises from the previous stabil- 
ity analysis, refers to the final fate of the unstable MSBS 
with rj > rjmax- We address this issue by performing 
long evolutions of unstable MSBS configurations, until 
they reach a stationary state. In order to accelerate the 
growth of the unstable modes, we perturbed the MSBS 
with a massless scalar field located far from Rg^. As 
explained in the previous section, this scalar field inter- 
acts gravitationally with the MSBS, perturbing it slightly 
and exciting the unstable modes. These modes grow ex- 
ponentially, starting with a small amplitude, result that 
can be obtained also from a linear perturbation analysis. 
When the amplitude of these perturbations is larger, the 
nonlinear effects become important and the evolution can 
only be followed numerically in order to discern the final 
state of the MSBS. 

Fig. m displays the number of particles in different 
states and the total number of particles for two MSBS 
with the same (/)i(0) = 0.0197 with a fraction given by 
77 = 3 and 77 = 0.4. The total Noether charge remains 
almost constant in the stable case 77 = 0.4, showing the 
accuracy of the numerical code within a 0.004% error 
in this quantity. The unstable case 7/ = 3 exhibits scalar 
field radiation during the evolution, producing a decrease 
of around 18% in the total Noether charge, as it can be 
seen in the convergence test presented in fig. [9] This ra- 
diation translates into a small change in the amplitude 
of the scalar fields. Taking a closer look at the maximum 
value of the scalar fields in the center and the frequen- 
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FIG. 6: The same as in figO for different resolutions and 
perturbation amplitudes. Despite a delay in the excitement 
of the exponential mode of the unstable case, the results are 
robust with respect to changes in resolution and perturbation. 



j(0) = 0.0143 
□ (|)j(0) = 0.0197 
(t),(0) = 0.0423 



2x10'^ 



IxlO 




FIG. 7: Fitting of the exponential growth for three families 
of MSBS configurations. The extrapolated value corresponds 
to the maximum allowed stable fraction. 



cies, displayed in figure [TOl one can notice a change in the 
position of the node; the excited state has decayed to a 
ground one, while the ground one has jumped to the first 
excited state. With this "flip-flop" of the scalar flelds, 
the final 77 is in the stable domain. A similar behavior 
is observed for all the unstable MSBS configurations in- 
cluded in the study, which indicates that this could be a 
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FIG. 8: The number of particles in the ground state and in 
the excited state for the fraction 77 = 3, together with the 
total number of particles for the fraction 77 = 3 and 77 = 0.4. 
There is a significant loss of the number of particles in the 
unstable configuration, while the stable only looses 0.004% 
due to numerical dissipation. 



common feature of their evolution. The time required by 
an unstable MSBS to settle down into a stable configu- 
ration, increases as the fraction gets closer to rjmax- 
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FIG. 9: The total number of particles for the fraction 7; = 3 
computed with three different resolutions. The loss in the 
number of particles converges to a value around 18%. 

The previous statements can be seen graphically in Fig. 
[TT] The figure shows the collection of initial configura- 
tions, each one characterized by the number of particles 
in the ground state A'^i , and the number of particles in 
the excited state N2- Every pair {Ni,N2) corresponds 
to a configuration with a pair of eigen-values {(^1,0:2), 
although only uji is shown in the figure. On the top of 
these collection of initial configurations, we have plotted 
the time evolution of the previous simulations. The red 
dot corresponds to the configuration labeled with fraction 
rj = 0.4 and the blue dots represent the time evolution for 
the configuration labeled with fraction 77 = 3.0. Apply- 
ing the same perturbation for the two different states, we 
can discern two different behaviors: the fraction 77 = 0.4 
is stable and remains static, while the configuration with 
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time 

FIG. 10: The maximum of tlie central value of tlie scalar 
fields in tlie different states for the fraction 77 = 3 on the top, 
and the frequencies of those modes on the bottom. There is 
a flip-flop of the scalar field and the frequencies a,t t = 5000. 

?; — 3.0 is unstable and evolves far from the original con- 
figuration (that lies on the sheet of static configurations) 
to a "forbidden" region. The jump is due to the flip of 
ground state -> excited state and viceversa. Then, it 
starts to lose scalar field (i.e., evolves in (7Vi,iV2)) and 
slowly approaches an equilibrium configuration on the 
sheet of equilibrium configurations. 

The final state of the MSBS with 77 = 3 can be inferred 
from the central value of the amplitude of the scalar fields 
(top panel Fig. I10[) . It is important to remind that now, 
due to the switch of the frequencies, the new ground state 
corresponds to (t>2{0) while the new excited state corre- 
sponds to 1^1(0). The amplitudes oscillate around a cen- 
tral value, namely 

0i(O)"^'^ = 02(0, i/w) = 0.058 ±0.004 
02(0)"^'" = 01 (0,t/w) =0.016 ±0.004. (34) 

Using these values as initial conditions for the new equi- 
librium configuration we can compute the eigenvalues 
and the number of particles shown in Table ID In the 
same table are displayed the values obtained from the 
evolution of the MSBS configuration with ij = 3. There 
is a remarkable agreement in the eigenvalues, but they 
still differ in the number of particles. This means that 
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Equilibrium 


0.87 


0.94 


0.604 


0.126 


Evolved A X = 0.020 


0.864 ± 0.004 


0.936 ±0.006 


0.71 


0.16 


Evolved A X = 0.015 


0.862 ± 0.003 


0.934 ±0.010 


0.706 


0.126 


Evolved A X = 0.010 


0.872 ±0.010 


0.940 ± 0.005 


0.696 


0.124 



TABLE I: Expected eigen-values and number of particles for 
the "new" equilibrium configurations. The values obtained 
from the late time evolution of the configuration with = 3 
is shown for comparison. 



the system will still loose particles at a slow rate during 
the evolution, as corroborated by Fig. [TTJ The slow loss 
of particles has been observed in finite perturbed system 
of excited states and even tough the final state has been 
inferred in the same way as we have done here (ssl . Issj . 
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FIG. 11: "Stable" and "unstable" states 



V. DARK MATTER HALOS 

In previous sections we have constructed MSBS with 
two states, a ground and a first excited state, and shown 
their stability. In this section, we come back to our initial 
motivation and we will illustrate how MSBS can lead to 
RCs which are in better agreement the RCs of galaxies 
in the context of SFDM. A detailed analysis of fitting the 
RC of galaxies including baryonic matter and experimen- 
tal data is out of the scope of the present work. Instead, 
we will just present a comparison between the behavior 
of a test particle immersed in the gravitational poten- 
tial produced by a MSBS and a single BS, showing that 
there is an improvement in the sense that the MSBS has a 
fiatter profile far from the center. Neglecting the baryonic 
contribution is a reasonable assumption in some galaxies 
such as the Low Surface Brightness Galaxies. 

For the static spherically symmetric metric considered 
here the circular orbit geodesic obeys (45| . 

z)^ = ra dra . (35) 
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FIG. 12: Density profiles and rotational curves for a MSBS 
with (^i(O) = 0.0197 and fraction 77 = 1, and the standard 
Boson-star made of a single state scalar field with the same 
central value <^i(0) = 0.0197 (77 = 0). Improvement in the 
keplerian tale is observed for large values of r. 



As an illustrative example, the top panel of Fig. [T^] 
presents a comparison between two rotational curves ob- 
tained for a single state boson star with (/)i(0) — 0.0197 
(ie, a MSBS with 77 = 0) and for a MSBS with the same 
amplitude 4>i{0) = 0.0197 of the ground state scalar field, 
and the same number of particles in the first excited state 
(ie, 77 = 1). The region with a flat plateau is larger, 
suggesting that for higher excited states (or MSBS with 
several higher states) the region with constant rotational 
velocities could be extended to larger radii. 

In order to understand this behavior, it will be helpful 
to see the mass density profile defined as 



p(r) 



1 dMjr) 

J.2 flj. 



(36) 



and it is shown in the bottom panel of Fig. [T^] for the 
same configurations mentioned above. We can see that 
for the single BS (77 = 0) the density decays exponentially 
as r — > c», making it difficult to fit the flat rotational 
curve profiles present in most galaxies. However, the 
MSBS configurations with large Noether fractions have 
a radius which is significantly larger than the one corre- 
sponding to the single ground state, with an exponential 
decay only in the tail of the excited state. 



Another issue related with boson stars as dark mat- 
ter models, was the lack of degrees of freedom to match 
the different sizes and masses of the observed galaxies. 
For a single boson star without self-interaction, the only 
free parameters are the mass of the boson particle m and 
the central value of the scalar field (/)(r — 0), which de- 
termines the compactness of the object (ie, ratio of total 
mass over radius) in adimensional units. There have been 
several attempts to fit these parameters [23-l28| with dif- 
ferent levels of success. By allowing more general MSBS, 
there are extra free parameters coming from the different 
fractions between the ground and excited states. These 
parameters change not only the total mass, but also the 
compactness of the final object. The extra degrees of 
freedom may allow a better fit of the models to different 
galaxies. 



VI. CONCLUSIONS 

We have constructed generalized boson star configu- 
rations, where two coexisting states of the scalar field 
are present. Our initial data construction is based on 
two main quantities that describe them: the gravita- 
tional mass M and the radius /igg. We have shown that 
these boson stars are stable under small radial perturba- 
tions, for a certain range of the fraction (77 < 1) between 
the Noether charges. These results may sound counter- 
intuitive, given the known fact that single BS in excited 
states are unstable under finite perturbations. Neverthe- 
less, the addition of an extra scalar field allows for an 
infinite number of new equilibrium configurations. The 
known plot of M vs. 0i(O) for ground state or excited 
state single BSs is now extended, getting a different curve 
for each fixed value of (^2(0) — constant^ as it is shown 
in the figure[T3l The single BS in the ground state corre- 
sponds to the case (^2(0) = 0, which only has an extreme 
at the maximum allowed mass. The configurations on the 
right of that point (marked with a triangle) are unstable. 
For low values of 02 (0) these curves contain now two ex- 
tremes. In addition to the maximum allowed mass, there 
is a minimum close to the fraction rj « 0.5. Although the 
presence of a new extreme in the curves could suggest 
a change in the stability around the fraction 77 « 0.5, 
our numerical stability analysis presented in the section 
IIV Al indicates that the dynamical stability regime is ex- 
tended beyond this point up to rjmax ~ f ■ Furthermore, 
we have found also stable and unstable configurations in 
regions with no extremes in the figure [T3l This supports 
the idea that simple stability arguments from the single 
boson star case are not valid in this case and that a more 
detailed future study is necessary in order to confirm and 
understand completely this problem. 

The unstable configurations evolve and settle down 
into stable configurations. MSBS allow to obtain a flat 
region in the velocity rotational curves, as shown by the 
examples in this paper. We considered cases with only 
two different states of the scalar field. As future work, we 
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FIG. 13: Mass as a function of (/f>i(0) for different values of 
<?f'2(0) = constant. The triangles mark the maximum allowed 
masses (ie, the maximum of the curves) while the circles corre- 
spond to the minimums of the curves. The maximum fraction 
Tjmax, displayed with squares, is found numerically on the left 
of the minimum of the curves. 



are planning to construct MSBS where several states are 
coexisting. These models allovif more degrees of freedom, 
and could be used to fit accurately the rotational curves 
within the observational data. 
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as a classical field, so the source in the r.h.s. of the 
Einstein equations has the expectation value of (IA2p over 
a state of the system of many particles jQ), namely 

Gab^&Tl{Q\fab\Q). (A3) 

The operators in the quantized field expansion of eq. 
(IA1[) can be interpreted as creation bnim and annihila- 
tion quantum operators. These operators satisfy 



the following commutation relations: 



I'm' 



0. 



(A4) 
(A5) 



The coefficients of the scalar field operator in eq. (|Aip 
must satisfy the Klein Gordon (KG) equation in a curved 
space time eq. ([5]). 

Using the relations (jA4p . one can construct the states 



\Q) 



\Nnl7n, Nn'l'm' .Nni'l'ir. 



(A6) 



These states are orthonormal and represent particle 
states, each composed of N scalar particles distributed 
in sets of Nnim particles of mass fj,, with angular momen- 
tum and azimuthal momentum hm^. The n subindex 
labels the energy eigenstate. Then the expectation value 
of the energy momentum operator in (jASp can be calcu- 
lated as 



{Tab) = {Nnlm,Nn'l' 



\Tab\-^nl'ni^ -^ri'l'm' i 



(A7) 

The orthonomafity of the quantum states ensures that 
this expectation value is given as a superposition of the 
expectation values of the energy-momentum for each 
state. Then (|A3P is given by 



Gab = ^T^'y^^Cnlm{Nnlm\Tab\Nnlm)^ (A8) 
nlm 

where Cnim are normalization coefficients [2^ . Therefore, 
in the case where more than one state is populated, the 
source of the Einstein equations is equivalent to the su- 
perposition of many uncoupled scalar fields. Each field 
generates its own stress energy tensor {Nnim\Tab\Nnim) ■ 



Appendix A: Equivalence of real quantized scalar 
field and multi boson stars 



The many boson-system is described by a second quan- 
tized free scalar field 

* = ^5„;™$„i„,(i,x) -I- 6t (t,x), (Al) 

nlm 

with an energy-momentum tensor operator given by 

Tab - da^db^ - ^^gabig^^OAdd^ + /i'l^n, (A2) 

where fx is the scalar field mass and the convention h — 
c = 1 has been adopted. The gravitational field is treated 



Appendix B: The Z3 system in spherical symmetry, 
normal coordinates and regularization 

The line element of a generic spherically symmetric 
spacetime can be written as 



ds^ = —a^dt^ + grrdr^ 



'geedfl^ 



(Bl) 



where we made explicit the singular factor r^, such that 
the metric components are regular. However, this change 
amounts to a transformation of the variables 



gee 



Dr- 



r'^gee, 
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where the quantities marked with tilde are the variables 
typically used in spherical symmetry. In order to ensure 
the stability of the implementation, one has to deal with 
the factors 1/r in the fluxes and 1/r^ in the sources. 



A regular system of evolution equations can be ob- 
tained, by ensuring a cross-cancellation between the sin- 
gular terms. We take advantage of the way the mo- 
mcntiim constraint was built into the system and rede- 
fine the variable Zr in order to obtain the desired cross- 
cancellation, 



4r 



1- 



9ee 



We can eliminate this way the singularities from the evo- 
lution variables and the numerical errors caused by the 
geometrical factors in the fluxes and sources. One can no- 
tice that the sources contain terms like 1/r times other 
variables which are radial derivatives of the metric coeffl- 
cients. But these terms do not create problems at r — )• 0, 
as the radial derivatives of any diffcrentiable function 
vanish at the origin. However, due to finite differencing, 
we can not use a grid point at r = 0. 



Klein-Gordon system in first order form is: 

dtSrr = —2agrrKr 



9 



dtAr = -dr[aftrK], 
dtDr/ = —dr[aK/], 
dtDre' = -dr[aKe\ 
dtZr = -dr[2aKe^] + 



- KJ 



Z. + J-fl-^ 
4r V gee 



4r gee 



dtK/ = —dr 



ag 



"^ 3 3 



-I- 



+ ^[g'-'^iDr/ -Ar- AZr) + g'^Dro' - Ar)] + 

ir 



3» 



ag 



+ 



+a 



^¥''{Ar - 2Dr/ - 4Zr) + g'^Ar - 2£»,/)]- 
4r 



{Dr/ - Dre^ - 2Ar)+ 
+ ^g'- (Ore' + {Dr/ - 4Ar)+ 



dt(l) = ay/g^(j)t, 

dt(i)r = dr[ayf^(l)t], 

dt4>t = dr[ay/^(t)r] + aVF^[2(£»^/ + l/r)(/>r + 
+2^/g/K0^ (j)t - m^grrcj)]- 

The complex scalar field can decomposed as 

= <^«-#^ 
6 = d^^ + id)' 



where cj)]^ is the real part, 0/ the imaginary part and (p 
The final set of equations for the regularized Einstein- its complex conjugate. 
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The matter terms can be explicitly written in terms of 
the components of the scalar field: 



So' = Wvtf+{^^f]-g'Vrr+{m- 



~M\^y + (0^)1} + -g^UtY - {AY 

where ip is the scalar field perturbation. 
The charge density can be computed as 



N = aJ° 



1 



-.{<l>'cf>f-cl>^cl>l). 



number of bosonic particles 

N = j Vh N dx^ = 4tt j N yj^geedr. 

We compute both the ADM and the Tolman masses in 
order to check our numerical evolutions. The ADM mass 
is defined as 

Madm = lim / gP^idggpk - dkgpq]N^dS, 

where A'''" = ^g"8/ is the unit outward normal to the 
sphere. In our coordinates, it can be translated into 



Madm = -r^VF'i^Sr/. 
The Tolman mass can be calculated as 

MtoI = JiTo" -T,')^ dx^ = 

= -4Trr^ay/g;;gee{T + S/ + 2Se'). 



(Bl) 



The space volume integral of N can be interpreted as the 
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